16. 混合ガウシアンモデル(Gaussian Mixture Model, GMM)

16.1. EMアルゴリズム

EMアルゴリズムは、EステップとMステップを交互に繰り返すことで、隠れ変数を含む確率モデルのパラメータを推定する手法である。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# ----------------------------
# データ生成 (Ground Truth)
# ----------------------------
np.random.seed(42)

true_means = np.array([[0, 0], [3, 3], [-3, 3]])
true_covs = [
    np.array([[1, 0.3], [0.3, 1]]),
    np.array([[1, -0.2], [-0.2, 1.5]]),
    np.array([[1.5, 0.4], [0.4, 1]])
]
true_weights = np.array([0.3, 0.4, 0.3])

n_samples = 10**3
X = np.vstack([
    np.random.multivariate_normal(mean, cov, int(w * n_samples))
    for mean, cov, w in zip(true_means, true_covs, true_weights)
])

# ----------------------------
# ガウス分布のPDF
# ----------------------------
def gaussian_pdf(x, mean, cov):
    d = x.shape[1]
    det = np.linalg.det(cov)
    inv = np.linalg.inv(cov)
    norm_const = 1.0 / np.sqrt((2 * np.pi) ** d * det)
    diff = x - mean
    return norm_const * np.exp(-0.5 * np.sum(diff @ inv * diff, axis=1))

# ----------------------------
# プロット: データと推定されたガウス分布の等高線
# ----------------------------
def plot_gaussian_ellipse(mean, cov, ax, color='red'):
    vals, vecs = np.linalg.eigh(cov)
    order = vals.argsort()[::-1]
    vals, vecs = vals[order], vecs[:, order]
    angle = np.degrees(np.arctan2(*vecs[:, 0][::-1]))
    width, height = 2 * np.sqrt(vals)
    ellip = Ellipse(xy=mean, width=width, height=height,
                    angle=angle, edgecolor=color, fc='None', lw=2)
    ax.add_artist(ellip)

def snapshot_gaussians(X, mu, sigma, K, iteration):
    plt.figure(figsize=(5, 5))
    plt.scatter(X[:, 0], X[:, 1], s=10, alpha=0.5)

    colors = ["red", "green", "blue"]
    for j in range(K):
        plot_gaussian_ellipse(mu[j], sigma[j], plt.gca(), color=colors[j])
        plt.scatter(mu[j, 0], mu[j, 1], c=colors[j], marker="x", s=100)    
    plt.title(f"Data and Fitted Gaussian Components @ Iteration = {iteration}")
    plt.xlabel("x1")
    plt.ylabel("x2")
    plt.show()


# ----------------------------
# EMアルゴリズム (GMM)
# ----------------------------
def EM_GMM(X, K, num_iters=100, tol=1e-5):
    n, d = X.shape
    np.random.seed(11)
    # 初期化
    phi = np.ones(K) / K
    mu = X[np.random.choice(n, K, replace=False)]
    sigma = np.array([np.cov(X, rowvar=False) for _ in range(K)])

    log_likelihoods = []
    snapshots = []

    for iteration in range(num_iters):
        # E-step
        responsibilities = np.zeros((n, K))
        for j in range(K):
            responsibilities[:, j] = phi[j] * gaussian_pdf(X, mu[j], sigma[j])
        responsibilities /= responsibilities.sum(axis=1, keepdims=True)

        # M-step
        N_k = responsibilities.sum(axis=0)
        phi = N_k / n
        mu = (responsibilities.T @ X) / N_k[:, np.newaxis]
        for j in range(K):
            diff = X - mu[j]
            sigma[j] = (responsibilities[:, j][:, np.newaxis] * diff).T @ diff / N_k[j]

        # 対数尤度
        ll = np.sum(np.log(np.sum([
            phi[j] * gaussian_pdf(X, mu[j], sigma[j]) for j in range(K)
        ], axis=0)))
        log_likelihoods.append(ll)

        # スナップショット保存
        snapshots.append((mu.copy(), sigma.copy(), iteration))

        # 収束判定
        if iteration > 0 and abs(log_likelihoods[-1] - log_likelihoods[-2]) < tol:
            print(f"Converged at iteration {iteration}")
            break

    # アニメーション作成

    fig, ax = plt.subplots(figsize=(5, 5))
    scatter = ax.scatter(X[:, 0], X[:, 1], s=10, alpha=0.5)
    title = ax.set_title("")

    def animate(i):
        ax.clear()
        ax.scatter(X[:, 0], X[:, 1], s=10, alpha=0.5)
        mu_snap, sigma_snap, iter_snap = snapshots[i]
        colors = ["red", "green", "blue"]
        for j in range(K):
            plot_gaussian_ellipse(mu_snap[j], sigma_snap[j], ax, color=colors[j])
            ax.scatter(mu_snap[j, 0], mu_snap[j, 1], c=colors[j], marker="x", s=100)
        ax.set_title(f"Data and Fitted Gaussian Components @ Iteration = {iter_snap}")
        ax.set_xlabel("x1")
        ax.set_ylabel("x2")

    anim = FuncAnimation(fig, animate, frames=len(snapshots), interval=500)
    plt.close(fig)
    display(HTML(anim.to_jshtml()))

    return phi, mu, sigma, responsibilities, log_likelihoods

# ----------------------------
# 実行
# ----------------------------
K = 3
phi, mu, sigma, responsibilities, log_likelihoods = EM_GMM(X, K)

print("Estimated mixing coefficients (phi):", phi)
print("Estimated means (mu):\n", mu)
print("Estimated covariances (sigma):\n", sigma)

# ----------------------------
# プロット: 対数尤度の推移
# ----------------------------
plt.figure(figsize=(6, 4))
plt.plot(log_likelihoods, marker="o")
plt.title("Log-Likelihood over EM Iterations")
plt.xlabel("Iteration")
plt.ylabel("Log-Likelihood")
plt.grid(True)
plt.show()

Converged at iteration 57

Once Loop Reflect

Estimated mixing coefficients (phi): [0.3103679  0.30013538 0.38949671]
Estimated means (mu):
 [[ 0.08871446  0.04564594]
 [-3.08442868  3.07309071]
 [ 3.05915799  3.15711097]]
Estimated covariances (sigma):
 [[[ 1.01691342  0.38486941]
  [ 0.38486941  0.98651172]]

 [[ 1.44451192  0.28724678]
  [ 0.28724678  0.91035083]]

 [[ 1.09349783 -0.24274844]
  [-0.24274844  1.27556857]]]

true_means = np.array([[0, 0], [3, 3], [-3, 3]])
true_covs = [
    np.array([[1, 0.3], [0.3, 1]]),
    np.array([[1, -0.2], [-0.2, 1.5]]),
    np.array([[1.5, 0.4], [0.4, 1]])
]
true_weights = np.array([0.3, 0.4, 0.3])